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Abstract 

We present an approach to solid-state electronic-structure calculations based 
on the finite-element method. In this method, the basis functions are strictly 
local, piecewise polynomials. Because the basis is composed of polynomials, 
the method is completely general and its convergence can be controlled sys- 
tematically. Because the basis functions are strictly local in real space, the 
method allows for variable resolution in real space; produces sparse, struc- 
tured matrices, enabling the effective use of iterative solution methods; and 
is well suited to parallel implementation. The method thus combines the 
significant advantages of both real-space-grid and basis-oriented approaches 
and so promises to be particularly well suited for large, accurate ab initio 
calculations. We develop the theory of our approach in detail, discuss advan- 
tages and disadvantages, and report initial results, including the first fully 
three-dimensional electronic band structures calculated by the method. 
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I. INTRODUCTION 



Over the course of the last few decades, the density functional theory (DFT) of Hohen- 
berg, Kohn, and ShamS has proven to be an accurate and reliable basis for the ab initio 
calculation of a wide variety of materials properties. But the solution of the Kohn-Sham 
equations is a formidable task, and this has significantly limited the range of physical systems 
which can be considered. 

Among the most popular methods of solving the equations has been the plane-wave 
(PW) method — typically coupled with pseudopotentials to eliminate core electrons.! For all 
its advantages, however, the PW method has some notable disadvantages with respect to 
the solution of large problems. First, because the PW basis functions are not local in real 
space, they give rise to a dense Hamiltonian matrix which in turn limits the effectiveness 
of iterative solution methods.! Second, the method requires Fourier transforms which are 
difficult to implement efficiently on massively parallel architectures due to the need for 
non-local communications. Finally, the PW basis has the same resolution at all points in 
space, which causes considerable difficulties in the treatment of highly localized systems 
such as first-row elements and transition metals. Recent progress on this problem has 
included ultrasoft pseudopotentials,@ optimized pseudopotentials,!'! and adaptive coordinate 
transformations.ETH 

The limitations of the PW approach have inspired the development of various 
"real-space" approaches including finite-difference (FD),lll~!3 finite-element (FE),ii~§l and 
wavelet&H methods. Of these, perhaps the most mature and successful to date have been 
the FD methods. These methods produce sparse, structured Hamiltonian (and in some 
cases, overlap) matrices, require no Fourier transforms, and allow for some degree of vari- 
able resolution in real space. But to achieve these advantages, these methods give up the 
use of a basis and work instead by discretizing individual terms of the differential equation 
of interest on a real-space grid. As a result, quantities of interest are defined only at a 
discrete set of points in space, limiting the accuracy of integrations and complicating the 
handling of singular functions such as all-electron potentials. And as a further consequence, 
the methods are not variational: the error can be of either sign and convergence is often 
from below. 

Finite-element methods&H achieve the significant advantages of FD methods without 
giving up the use of a basis. Like the PW method, the FE method is an expansion method. 
In the FE method, however, the basis functions are strictly local, piecewise polynomials. A 
simple 1-D example is shown in Fig. [I] (which we discuss further below). Because the basis is 
composed of polynomials, it is completely general and the convergence of the method can be 
controlled systematically by increasing the number or order of basis functions. Because the 
basis functions are strictly local in real space, the method achieves the significant advantages 
of FD approaches: The method produces sparse, structured matrices, which in turn enable 
the effective use of iterative solution methods. The method requires no Fourier transforms, as 
all calculations are performed in real space. And the method allows for variable resolution 
in real space — more so than FD approaches — by increasing the number or order of basis 
functions where needed. The method thus combines the significant advantages of both 
real-space-grid and basis-oriented approaches. 

Some disadvantages of FE methods are that the matrices produced tend to be less sparse, 
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and often less simply structured, than those produced by FD methods, and that these 
matrices must be stored. In addition, FE methods produce generalized rather than standard 
eigenvalue problems, as produced by many FD approaches, and they can be significantly 
more difficult to implement than FD or PW approaches. 

The FE methodli^rEil has had a long history of success in quite diverse applications ranging 
from civil engineering to quantum mechanics. Applications in engineering go back to the 
1950s. Applications to the electronic-structure of isolated atomic and molecular systems 
began to appear as early as the 1970s£3 White, Wilkins, and Tetei§ applied the method 
to full 3-D atomic and molecular calculations in 1989 and demonstrated the advantageous 
scaling of the method with the number of basis functions, afforded by the strict locality and 
real-space nature of the basis. 

There have, however, been relatively few applications to solids. Hermansson and YevickEl 
applied the FE method to full 3-D solid-state electronic-structure calculations in 1986, but 
having reached a negative conclusion in the study of small systems with uniform meshes, 
discussed their approach only generally and, though apparently capable of arbitrary Bril- 
louin zone sampling, reported only T-point results. More recently, Tsuchida and Tsukada@ 
have applied the method to full 3-D molecular and solid-state electronic-structure calcula- 
tions. They have implemented self-consistency, non-uniform meshes, adaptive coordinate 
transformations, pseudopotential and all-electron calculations, and have demonstrated the 
favorable efficiency of the method relative to FD approaches. Their solid-state results have, 
however, also been limited to the T-point. 

Here, we present a full 3-D FE approach to solid-state electronic-structure calculations 
which allows arbitrary sampling of the Brillouin zone and report, to our knowledge, the first 
such band structures calculated by the method. Our development differs from that of Ref. 25 
in that we have taken a Galerkin approach.0 Our development is closer to that of FerrarijH 
who developed a Galerkin approach allowing arbitrary sampling of the Brillouin zone in the 
context of 2-D super-lattice calculations. As we have not yet implemented self-consistency, 
our results here are limited to model potentials and empirical pseudopotentials. 

The remainder of the paper is organized as follows: In Sec. |T| we discuss the details 
of our approach. We begin in Sec. [11 A| with a description of the basis. In Sec. [11 B| we 
show how the basis can be applied to the solution of the Schrodinger equation subject to 
boundary conditions appropriate to a periodic solid. In Sec. JT| we present results for a model 
potential and Si pseudopotential, including band structures and details of the convergence 
of the method. The conclusions in Sec. [TV] are followed by an appendix giving the details of 
the particular 3-D basis which we employ. 



II. METHOD 



A. Basis 



Finite-element bases consist of strictly local, piecewise polynomials. They are constructed 
generally as follows: The domain is partitioned into subdomains called elements. Within 
each element a set of polynomial basis functions is defined. These element polynomials are 
then pieced together at inter-element boundaries to form the piecewise polynomial basis 
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functions of the method. In order to apply the method to periodic problems, we take the 
additional step here of piecing together element polynomials across the domain boundary.il 

The essential ideas are perhaps best conveyed by a simple example: a 1-D, periodic, 
piecewise-linear basis. Figure 0(a) shows the complete basis and Figs. [j](b)-(e) show the 
individual basis functions. In this case, the domain [a, b] has been partitioned into four 
elements. For simplicity we have defined a uniform partition but this need not be the case 
in general. Within each element, two linear_polynomial basis functions have been defined 
to make the basis complete to linear order.a More and higher-order polynomials can be 
used to increase the order of completeness. Figure |I](b) shows the basis function which 
results from piecing together element polynomials across the domain boundary. Figures 
|l](c)-(e) show the basis functions which result from piecing together element polynomials 
across inter-element boundaries. 

Note that the resulting basis functions are C°-continuous, i.e., continuous but not nec- 
essarily smooth.il Smoother bases can be constructed by requiring continuity of higher 
derivatives (which would require higher-order polynomials), but a C° basis offers a unique 
and potentially significant advantage: an efficient and natural representation of the wave- 
function and charge density cusps which occur in all-electron calculations and which cause 
such difficulty for conventional, necessarily smooth bases — and for FD approaches as well, 
as they also assume some degree of derivative continuity. We have therefore chosen a C° 
basis for our approach. 

The application of such a basis, however, to the solution of a second-order differen- 
tial equation, with periodic boundary conditions in particular, requires some consideration. 
First, the application of the Laplacian to such functions is clearly problematic. Second, 
referring again to Fig. |], note that the basis is value-periodic, i.e., 

<Ma + )=<M&-), (1) 

but not derivative-periodic, i.e., it is not the case for all basis functions that 

<P' l {a+)=<P> l {b-). (2) 

Thus, the satisfaction of periodic boundary conditions is nontrivial. We address these issues 
in Sec. Ill B, 



Referring again to Fig. [l|, note also that the basis functions take on a value of one at 
their associated nodes and zero at all others, i.e., 

4>i(xj) = 5ij. (3) 

Thus, a FE expansion, 

f{x)=^2ci^i{x), (4) 

is such that 

f(xi) = (H, (5) 
giving the expansion coefficients a direct real-space meaning. 
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Finally, and perhaps most significantly, we note that the basis functions are strictly local 
in real space, i.e., nonzero over only a (typically small) fraction of the domain. It is this 
property of the basis which allows the method to achieve the significant advantages of FD 
approaches. 

In our calculations, we have employed a 3-D, C°, piecewise-cubic basis. Many other 
choices are possible. Higher-order completeness generally leads to smaller matrices and 
higher-order convergence, but also to less sparseness. The details of the basis are given in 
the Appendix. 

B. Discretization 

We solve 

- V 2 ^ + Vip - sip = (6) 

in a unit cell, where V is an arbitrary periodic potential, as appropriate for a periodic solid. 

We begin by reducing the Bloch-periodic problem to a periodic one. Since V is periodic, 
we can take ip to be of the form 

^(x) = n(x)e ik - x , (7) 
where u is a complex, cell-periodic function satisfying 

m(x) = u(x + R) (8) 
for all lattice vectors R. Substitution of the form ([?]) into (0) then gives 

-X7 2 u-2ik-Vu + (V + k 2 -e)u = 0. (9) 
From the periodicity condition (^), we take the boundary conditions to be: 

w(x) = w(x + Ri) Vx G Ti, 1 = 1,2, 3, (10) 

and 



n • Vu(x) = n ■ Vu(x + R z ) Vx G r,, 1 = 1, 2, 3, 



(11) 



where Ti and R; are the surfaces of the boundary T and associated lattice vectors shown in 
Fig. |2], and n is the outward unit normal at x. We denote the domain by Q and take it to 
be a parallelepiped for definiteness. The problem is thus reduced to Eqs. (P)- (|TTD . 

To facilitate the use of a C° basis, we next derive an equivalent "variational formulation" 
of the problem. The inner product of the differential equation (|9|) with an arbitrary "test 
function" v gives 



-V 2 m - 2i k ■ V« + (V + k 2 - e)u 



dQ = 0. 



(12) 
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Since v is arbitrary, the integral equation ([TJ) is equivalent to the differential equation @. 
To reduce the order of the highest derivative and produce a boundary term (whose usefulness 
will become clear subsequently), we integrate the V 2 term by parts:ll 



v*Vu ■ ndT 



-2i k • Vu + {V + k 2 - e)u 



dn = o. 



(13) 



To incorporate the "natural" boundary condition (PD, we now restrict v to 

v e V = {v : u(x) = u(x + Rj) Vx G Z = 1, 2, 3} , (14) 

i.e., to satisfy the "essential" boundary condition fllQ|) . Then, using the fact that the domain 
is a parallelepiped, the boundary term can be written as 



J v*(x) [V«(x) - Vm(x + R,)] • hdT, 



which vanishes upon the assertion of the natural boundary condition ([11]) . Thus, with the 
restriction (P^), the differential equation and natural boundary condition together imply the 
integral equation, 



Vv* -Vudtt + 



-2i k • Vu + (V + A; 2 - e)u 



dVL = 0. 



(15) 



n 



Finally, using again the arbitrariness of v, it can be shown that the converse also holds, and 
thus that the differential formulation (|9|)-([TT|) is in fact equivalent to the following variational 
formulation: Find the scalars e and functions u G V such that 



Vv* -VudQ + 



-2i k ■ Vu + (V + k 2 - e)u 



dfl = \/v G V. 



(16) 



We have thus reformulated the original problem in such a way that (1) the highest derivative 
which occurs is of order one, and (2) only the essential boundary condition remains — the 
natural boundary condition having been built into the equation itself. The problem is thus 
now in a form which is suitable for approximate solution in a C° FE basis since all terms are 
well defined for such functions and since such a basis can be readily constructed to satisfy 
the required value-periodicity (e.g. Fig. []]). 

To find an approximate solution, we now restrict v and u to a finite-dimensional subspace 
V n C V. The problem is then reduced to: Find the scalars e and functions u G V n such 
that 



Vv* -VudVL + 



-2% k • Vu + {V + k 2 - e)u dVt = Vt> G V n . 



(17) 



We proceed by constructing a real C° FE basis, <pi ■ ■ ■ <f) n , which satisfies the remaining 
essential boundary condition and so spans a subspace V n C V (e.g. Fig. [l]). We then 
express u as a complex linear combination, 
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u = 2 c i^' 



(18) 



so that « G V n . Substitution of the expansion (|18|) , and the fact that (|T7|) is satisfied 
Vf G V n if it is satisfied for v — </>j, i = 1 . . . n, leads finally to a generalized eigenproblem 
for the coefficients Cj and eigenvalues e determining the approximate eigenf unctions and 
eigenvalues of the variational formulation, and thus of the original problem: 



He = eSc, 



(19) 



where 



H ij = J [V0< ■ - 2i k • 0jV0j + (y + A; 2 ) 0^- 



(20) 



and 



5*; 



(21) 



As in the PW method, given the expansion of the potential, the above matrix elements can 
be evaluated exactly, due to the polynomial nature of the basis. As in the FD method, the 
above matrices are sparse and structured, due to the strict locality of the basis. 



III. RESULTS 

We have tested our approach in a number of applications ranging from the band structure 
of Si to positron charge distributions in Ceo- Here we present results which demonstrate the 
accuracy and convergence of the method for a model potential, for which analytic results 
are available, and for the more physically interesting case of Si. 

Figures |3] and [| show results for a 3-D generalized Kronig-Penney model potential: 

V = Vid(x) + V lD (y) + V 1D (z), (22) 

where 

Vid(0 = (periodically repeated). (23) 

Figure [3] shows the band structure obtained with a 6x6x6 uniform mesh vs. analytic results at 
selected k-points, for Vo = 6.5 Ry, a = 2 a.u., and 6 = 3 a.u. More quantitative information is 
displayed in Fig. [| which shows the convergence of the fractional error {E FEM —E exact ) / E exact 
of the first few eigenvalues with increasing numbers of elements, at an arbitrary k-point. The 
variational nature of the method is clearly demonstrated: the errors are strictly positive and 
monotonically decreasing. The consistent, sextic convergence of the method is also clearly 
demonstrated: the asymptotic slope of ~ —6 on the log-log scale corresponds to an error of 
0(h 6 ), where h is the mesh spacing, consistent with FE asymptotic convergence theorems 
for the cubic-complete case.0 
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Figures [5] and |] show results for a Si pseudopotential.Ea Since our approach allows for 
the direct treatment of an arbitrary parallelepiped domain, we show results for a two-atom 
primitive cell. In contrast, recent FD approaches have been limited to a small subset of 
Bravais lattices and have reported only supercell results for Si. Figure |5| shows the sequence 
of band structures obtained for 3x3x3, 4x4x4, and 6x6x6 uniform meshes vs. exact values 
at selected k-points. (Here "exact values" are from a highly converged PW calculation, 
using a 54 Ry cutoff.) The variational nature and rapid convergence of the method are 
again clearly demonstrated. Also apparent in the very coarse 3x3x3 results, are the inexact 
degeneracies of certain eigenvalues at high-symmetry k-points: for example, the splitting at 
the T-point of the triply degenerate value at the top of the valence band and the splitting 
at the X-point of the doubly degenerate lowest value. As noted in Ref. |24|, this is due to the 
fact that the basis is not constrained to have the full symmetry of the crystal. Thus, to the 
extent that the eigenvalues are approximate, so are the degeneracies; and as the eigenvalues 
converge, the degeneracies become exact. By the 6x6x6 mesh, the splittings are no longer 
apparent. Figure |6| shows the convergence of the fractional error of the first few eigenvalues 
with increasing numbers of elements, at an arbitrary k-point. The variational nature and 
consistent, sextic convergence of the method are again clearly demonstrated. 



IV. SUMMARY AND CONCLUSIONS 

We have presented an approach to solid-state electronic-structure calculations based on 
the finite-element method. In Sec. Ill A| we discussed the details of the basis, the most 



important being its polynomial composition and strict locality, leading to its generality and 
suitability for large-scale calculations. In Sec. [II B| we developed our approach to the solution 
of the Schrodinger equation, subject to boundary conditions appropriate to a periodic solid, 
using a C° finite-element basis: yielding a general method for solid-state electronic-structure 
calculations, allowing arbitrary sampling of the Brillouin zone. In Sec. [TTT] we presented initial 
results illustrating the accuracy and convergence characteristics of the method in electronic 
band-structure calculations. The consistent, sextic convergence and variational nature of 
the method were demonstrated. 

The finite-element method combines the significant advantages of both real-space-grid 
and basis-oriented approaches and so promises to be particularly well suited for large, ac- 
curate ab initio calculations. The results to date are promising, but the application of the 
finite-element method to solid-state electronic-structure calculations is still in its infancy 
and whether it will ultimately prove superior to other approaches will only be known after 
much further development. Our approach has already proven effective in large (863 atoms), 
non-self-consistent positron distribution and lifetime calculations.!^ Work on the addition of 
self-consistency, optimization of numerical methods, and parallelization is underway. 
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APPENDIX: 

In this appendix we discuss the details of the 3-D FE basis used in this work. We have 
employed standard 3-D 32-node "serendipity" elements@ These afford C° flexibility and 
cubic completeness with a minimum number of basis functions per element. As in standard 
FE references, we list below the parent basis functions defined on the parent element: [—1, l] 3 . 
The basis functions associated with any particular element are derived from these by a 
transformation.0 This permits the construction of quite general element meshes, permitting 
the precise concentration of degrees of freedom in real space where needed. For simplicity, 
we have limited our implementation to affine transformations. These are general enough to 
permit the direct treatment of an arbitrary parallelepiped domain and thus of any Bravais 
lattice. 

The parent element and associated nodal positions (denoted by open circles) are shown 
in Fig. 0. The 32 parent basis functions 0, and associated nodes (Ci,ViXi) are listed below: 

(&» Vi, d) <k 

(±1, ±1, ±1) J(T + £ )(1 + r/ )(l + Co){9(f + v 2 + C) ~ 19} 
(±±, ±1, ±1) Ml - £ 2 )(1 + 9&)(1 + ^o)(l + Co) 
(±1, ±|, ±1) J(l - 77 2 )(1 + 9%)(1 + £o)(l + Co) 
(±1, ±1, ±|) |(1 - C 2 )(l + 9Co)(l + £o)(l + vo) 

where 

£o = &f , Vo = ViV, Co = CiC- 

Thus, for example, the parent basis function associated with the node (— ~, 1, 1) is 

^(1 — £ 2 )(1 — 3£)(1 + r])(l + C); the one associated with the node (—1, |, 1) is 

^(1 — r/ 2 )(l + 37/) (1 — £)(1 + C), etc. Each takes on a value of one at its associated node 

and zero at all others. 

Upon piecing together element basis functions across inter-element and domain bound- 
aries, the resulting periodic piecewise-cubic basis contains seven basis functions per element. 



9 



REFERENCES 



1 P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964); W. Kohn and L.J. Sham, Phys 
Rev. 140, A1133 (1965). 

2 See e.g., W.E. Pickett, Comp. Phys. Rep. 9, 115 (1989) for a comprehensive review. 

3 Matrix- vector multiplies — the central operations in iterative solution methods — scale at 
best as NlogN, as compared to N for sparse matrices, where N is the dimension of the 
matrices. 

4 D. Vanderbilt, Phys. Rev. B 41, 7892 (1990); K. Laasonen, R. Car, C. Lee, D. Vanderbilt, 
Phys. Rev. B 43, 6796 (1991). 

5 A.M. Rappe, K.M. Rabe, E. Kaxiras, and J.D. Joannopoulos, Phys. Rev. B 41, 1227 
(1990). 

6 J.S. Lin, A. Qteish, M.C. Payne, and V. Heine, Phys. Rev. B 47, 4174 (1993). 

7 F. Gygi, Europhys. Lett. 19, 617 (1992); Phys. Rev. B 48, 11692 (1993); 51, 11190 (1995). 

8 A. Devenyi, K. Cho, T.A. Arias, and J.D. Joannopoulos, Phys. Rev. B 49, 13373 (1994). 
9 D.R. Hamann, Phys. Rev. B 51, 7337 (1995); 51, 9508 (1995); 56, 14979 (1997). 

°J.R. Chelikowsky, N. Troullier, and Y. Saad, Phys. Rev. Lett. 72, 1240 (1994); J.R. 

Chelikowsky, N. Troullier, K. Wu, and Y. Saad, Phys. Rev. B 50, 11355 (1994); J.R. 

Chelikowsky, Phys. Rev. B 57, 3333 (1998). 
1 E.L. Briggs, D.J. Sullivan, and J. Bernholc, Phys. Rev. B 52, R5471 (1995); 54, 14362 

(1996); J. Bernholc, E.L. Briggs, D.J. Sullivan, C.J. Brabec, M. Buongiorno Nardelli, K. 

Rapcewicz, C. Roland, M. Wensell, Int. J. Quant. Chem. 65, 531 (1997). 
2 G. Zumbach, N.A. Modine, and E. Kaxiras, Solid State Comm. 99, 57 (1996); N.A. 

Modine, G. Zumbach, and E. Kaxiras, Phys Rev. B 55, 10289 (1997). 
3 F. Gygi and G. Galli, Phys. Rev. B 52, R2229 (1995). 

4 K.A. Iyer, M.P Merrick, and T.L. Beck, J. Chem. Phys. 103, 227 (1995); T.L. Beck, K.A. 

Iyer, and M.P. Merrick, Int. J. Quant. Chem. 61, 341 (1997); T.L. Beck, Int. J. Quant. 

Chem. 65, 477 (1997). 
5 T. Hoshi, M. Arai, and T. Fujiwara, Phys. Rev. B 52, R5459 (1995). 
6 D.S. Burnett, Finite Element Analysis (Addison- Wesley, Reading, 1987). 
7 O.C. Zienkiewicz and R.L. Taylor, The Finite Element Method, 4th. ed. (McGraw-Hill, 

London, New York, 1988). 
8 K.-J. Bathe, Finite Element Procedures (Prentice Hall, Englewood Cliffs, 1996). 

9 G. Strang and G.J. Fix, An Analysis of the Finite Element Method (Prentice-Hall, Engle- 
wood Cliffs, 1973). 

;0 K. Eriksson, D. Estep, P. Hansbo, and C. Johnson, Computational Differential Equations 

(Cambridge University Press, Cambridge, 1996). 
;1 B.D. Reddy, Introductory Functional Analysis ( Springer- Verlag, New York, 1998). 
!2 For a review, see J. Linderberg, Comp. Phys. Rep. 6, 209 (1987); J. J.S. Neto and J. 

Linderberg, Comp. Phys. Comm. 66, 55 (1991). 
;3 S.R. White, JW. Wilkins, and M.P. Teter, Phys. Rev. B 39, 5819 (1989). 
;4 B. Hermansson and D. Yevick, Phys. Rev. B 33, 7241 (1986). 

;5 E. Tsuchida and M. Tsukada, Solid State Comm. 94, 5 (1995); Phys. Rev. B 52, 5573 

(1995); Phys. Rev. B 54, 7602 (1996). 
;6 J.E. Pask, B.M. Klein, and CY. Fong, Bull. Am. Phys. Soc. 43, 40 (1998). 



10 



27 



See e.g., Ref. 16 



28 R.L. Ferrari, Int. J. Numerical Modelling: Electronic Networks, Devices and Fields 6, 283 
(1993). 

29 K. Cho, T.A. Arias, J.D. Joannopoulos, and P.K. Lam, Phys. Rev. Lett. 71, 1808 (1993). 

30 S. Wei and M.Y. Chou, Phys. Rev. Lett. 76, 2650 (1996). 

31 C.J. Tymczak and X. Wang, Phys. Rev. Lett. 78, 3654 (1997). 

32 R.A. Lippert, T.A. Arias, and A. Edelman, J. Comp. Phys. 140, 278 (1998). 

33 For details see e.g., Ref. ch. 11; Ref. [16], ch. 5 and 13. The details of the additional 
step of piecing together across the domain boundary then follow straightforwardly from 
Sec. 0. 

34 A polynomial basis which is complete to order n spans the space of polynomials of order 
n. 

35 A function is of class C n if the function and its first n derivatives are continuous. 

36 We assume that v is sufficiently well-behaved to keep the integrals well-defined. For details 



regarding the relevant function spaces, see e.g., Refs. [19] and ^1 
37 See e.g., Ref. pi p. 426. 



38 M.L. Cohen and T.K. Bergstresser, Phys. Rev. 141, 789 (1966). 

39 P.A. Sterne, J.E. Pask, and B.M. Klein, Lawrence Livermore National Laboratory Report 
No. UCRL JC-131338 (to be published in Appl. Surf. Sci.) 

40 See e.g., Finite Element Handbook, edited by H. Kardestuncer, D.H. Norrie, et al. 
(McGraw-Hill, New York, 1987), p. 2.126. 

41 This is among the most common methods of generating FE bases. For details, see e.g., 



Ref. IB, ch. 8 and 13. 



11 



FIGURES 



(a) 



(b) 



(c) 



(d) 



(e) 




FIG. 1. Simple periodic finite-element basis: 1-D, piecewise-linear ) Basis, (b)-(e) 

Individual basis functions. The domain [a, b] is partitioned into elements (subdomains) (l)-(4) 
within which the basis functions are simple linear polynomials. The basis is thus simultaneously 
polynomial and strictly local in nature. 



R; 



FIG. 2. Parallelepiped unit cell (domain) $7, boundary T, surfaces Fi-t^, and associated 
lattice vectors R4-R3. 
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FIG. 3. FEM and exact band structures for 3-D generalized Kronig-Penney potential. FEM 
results are for a 6x6x6 uniform mesh of C° cubic elements. Exact results are from an analytic 
solution. 
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FIG. 4. Convergence of first few FEM eigenvalues for 3-D generalized Kronig-Penney poten- 
tial with increasing numbers of elements, at an arbitrary k-point. The convergence from above 
demonstrates the variational nature of the method. The asymptotic slope of ~ —6 demonstrates 
the sextic convergence of the method. 
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FIG. 5. FEM and exact band structures for Si pseudopotential. FEM results are for 3x3x3, 
4x4x4, and 6x6x6 uniform meshes of C° cubic elements. Exact results are from a highly converged 
plane-wave calculation. The rapid convergence and variational nature of the method are again 
demonstrated, with excellent agreement for the 6x6x6 mesh. 
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Iog 10 ( Elements in each direction) 
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FIG. 6. Convergence of first few FEM eigenvalues for Si pseudopotential with increasing num- 
bers of elements, at an arbitrary k-point. The variational nature and consistent, sextic convergence 
of the method are again demonstrated. 
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FIG. 7. Three-dimensional C° cubic parent element and associated nodal positions (denoted 
by open circles). 
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